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Abstract. In this paper we introduce Enzo, a 3D MPI-parallel Eulerian block- 
structured adaptive mesh refinement cosmology code. Enzo is designed to sim- 
ulate cosmological structure formation, but can also be used to simulate a wide 
range of astrophysical situations. Enzo solves dark matter N-body dynamics 
using the particle-mesh technique. The Poisson equation is solved using a 
combination of fast fourier transform (on a periodic root grid) and multigrid 
techniques (on non-periodic subgrids). Euler's equations of hydrodynamics 
are solved using a modified version of the piecewise parabolic method. Several 
additional physics packages are implemented in the code, including several 
varieties of radiative cooling, a metagalactic ultraviolet background, and pre- 
scriptions for star formation and feedback. We also show results illustrating 
properties of the adaptive mesh portion of the code. Information on profiling 
and optimizing the performance of the code can be found in the contribution 
by James Bordner in this volume. 



1 Introduction 



In astrophysics in general, and cosmology in particular, any given object of 
interest can have many important length and time scales. An excellent exam- 
ple of this is the process of galaxy formation. When studying the assembly of 
galaxies in a cosmological context, one wants to resolve a large enough volume 
of the universe to capture enough large-scale structure (a box with length on 
the order of several megaparsecs 5 ). However, in order to adequately resolve 

5 1 parsec = 3.26 light years = 3.0857 x 10 18 cm 
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structure in an individual galaxy one wants to have resolution two orders of 
magnitude smaller than the ultimate size of the objects of interest (a dwarf 
galaxy is on the order of ~ 1 kiloparsec). This typical cosmological problem 
requires roughly five orders of magnitude of dynamical range, which is pro- 
hibitively expensive when done using a single grid. Many other astrophysical 
phenomena, such as the study of molecular clouds and the formation of stars 
and galaxy clusters, require similarly large dynamical range. Many scientists 
have adopted Lagrangean techniques such as smoothed particle hydrodynam- 
ics (SPH)Q] to address these issues. However, this type of method suffers from 
several drawbacks, including poor shock resolution and fixed mass resolution 
in regions of interest. The use of grid-based techniques with structured adap- 
tive mesh refinement avoids many of these problems, and additionally allows 
the use of higher-order hydrodynamics schemes. 

In this paper we present Enzo, an MPI-parallel 3D Eulerian adaptive mesh 
refinement code. Though it was originally designed to study cosmological 
structure formation, the code is extensible and can be used for a wide range 
of astrophysical phenomena. For more information on the performance and 
optimization of the Enzo code, see the contribution by James Bordner in this 
volume. The Enzo web page, which contains documentation and the source 
code, can be found at http://cosmos.ucsd.edu/enzo/ 



2 Methodology 

2.1 Application and AMR Implementation 

Enzo is developed and maintained by the Laboratory for Computational As- 
trophysics at the University of California in San Diego. The code is written in 
a mixture of CH — h and Fortran 77. High-level functions and data structures are 
implemented in CH — h and computationally intensive lower-level functions are 
implemented in Fortran. Enzo is parallelized using the MPI message-passing 
library 6 and uses the HDF5 7 data format to write out data and restart files 
in a platform-independent format. The code is quite portable and has been 
ported to numerous parallel shared and distributed memory systems, includ- 
ing the IBM SPs and p690 systems, SGI Origin 2000s and numerous Linux 
Beowulf-style clusters. 

The code allows hydrodynamic and N-body simulations in 1, 2 and 

3 dimensions using the structured adaptive mesh refinement of Berger & 
ColellaQ, and allows arbitrary integer ratios of parent and child grid res- 
olution and mesh refinement based on a variety of criteria, including baryon 
and dark matter overdensity or slope, the existence of shocks, Jeans length, 
and cell cooling time. The code can also have fixed static nested subgrids, 

6 http : //www-unix.mcs . anl . gov/mpi/ 

7 |http : / /hdf . ncsa. uiuc . edu/HDF5/ 
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allowing higher initial resolution in a subvolume of the simulation. Refine- 
ment can occur anywhere within the simulation volume or in a user-specified 
subvolume. 

The AMR grid patches are the primary data structure in Enzo. Each in- 
dividual patch is treated as an individual object, and can contain both field 
variables and particle data. Individual patches are organized into a dynamic 
distributed AMR mesh hierarchy using arrays of linked lists to pointers to 
grid objects. The code uses a simple dynamic load-balancing scheme to dis- 
tribute the workload within each level of the AMR hierarchy evenly across all 
processors. 

Although each processor stores the entire distributed AMR hierarchy, not 
all processors contain all grid data. A grid is a real grid on a particular pro- 
cessor if its data is allocated to that processor, and a ghost grid if its data 
is allocated on a different processor. Each grid is a real grid on exactly one 
processor, and a ghost grid on all others. When communication is necessary, 
MPI is used to transfer the mesh or particle data between processors. The 
tree structure of a small illustrative 2D AMR hierachy - six total grids in a 
three level hierarchy distributed across two processors - is shown on the left 
in Figure H 



Distributed hierarchy 



Grid zones 
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Fig. 1. Real and ghost grids in a hierarchy; real and ghost zones in a grid. 



Each data field on a real grid is an array of zones with dimensionality 
equal to that of the simulation (typically 3D in cosmological structure forma- 
tion) . Zones are partitioned into a core block of real zones and a surrounding 
layer of ghost zones. Real zones are used to store the data field values, and 
ghost zones are used to temporarily store neighboring grid values when re- 
quired for updating real zones. The ghost zone layer is three zones deep in 
order to accomodate the computational stencil in the hydrodynamics solver 
fSection l2.3|) . as indicated in the right panel in Figure ^ These ghost zones 
can lead to significant computational and storage overhead, especially for the 
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smaller grid patches that are typically found in the deeper levels of an AMR 
grid hierarchy. 

For more information on Enzo implementation and data structures, see 
references 0, 0, and [|. 

2.2 N-body Dynamics 

The dynamics of large-scale structures are dominated by "dark matter," which 
accounts for ~ 85% of the matter in the universe but can only influence 
baryons via gravitational interaction. There are many other astrophysical sit- 
uations where gravitational physics is important as well, such as galaxy col- 
lisions, where the stars in the two galaxies tend to interact in a collisionlcss 
way. 

There are multiple ways that one can go about calculating the gravita- 
tional potential (which is an elliptical equation in the Newtonian limit) in 
a structured AMR framework. One way would be to model the dark matter 
(or other collisionless particle-like objects, such as stars) as a second fluid in 
addition to the baryon fluid and solve the collisionless Boltzmann equation, 
which follows the evolution of the fluid density in both physical space and 
velocity space (referred to collectively as "phase space" . This is computation- 
ally prohibitive due to the large dimensionality of the problem and because 
the interesting portion of the solution to the equation does not tend to oc- 
cupy a small volume of the computational domain, which makes this approach 
unappealing in the context of an AMR code. 

Enzo uses a totally different approach to collisionless systems, namely, the 
N-body method. This method follows trajectories of a representative sample 
of individual particles and is much more efficient than a direct solution of the 
Boltzmann equation in most astrophysical situations. The particle trajectories 
are controlled by a simple set of coupled equations (for simplicity, we omit 
cosmological terms) : 

If— « 

^ = -v* m 

Where x p and v p are the particle position and velocity vectors, respec- 
tively, and the term on the right-hand side of the second equation is the 
gravitational force term. The solution to this can be found by solving the 
elliptic Poisson's equation: 

V 2 = AnGp (3) 

where p is the density of both the collisional fluid (baryon gas) and the 
collisionless fluid (particles). 
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These equations are finite-differenced and for simplicity are solved with the 
same timestep as the equations of hydrodynamics. The dark matter particles 
are sampled onto the grids using the triangular-shaped cloud (TSC) interpo- 
lation technique to form a spatially discretized density field (analogous to the 
baryon densities used to calculate the equations of hydrodynamics) and the 
elliptical equation is solved using FFTs on the triply periodic root grid and 
multigrid relaxation on the subgrids. Once the forces have been computed on 
the mesh, they are interpolated to the particle positions where they are used 
to update their velocities. 

2.3 Hydrodynamics 

The primary hydrodynamic method used in Enzo is based on the piecewise 
parabolic method (PPM) of Woodward & Colella which has been signif- 
icantly modified for the study of cosmology. The modifications and several 
tests are described in Bryan et al. || , but we provide a short description here. 

PPM is a higher-order-accurate version of Godunov's method with third- 
order-accurate piecewise parabolic monotolic interpolation and a nonlinear 
Riemann solver for shock capturing. It does an excellent job capturing strong 
shocks and outflows. Multidimensional schemes are built up by directional 
splitting, and produce a method that is formally second-order-accurate in 
space and time and explicitly conserves energy, momentum and mass flux. The 
conservation laws for fluid mass, momentum and energy density are written 
in comoving coordinates for a Friedman-Robertson- Walker spacetime. Both 
the conservation laws and Riemann solver are modified to include gravity, as 
calculated in Section 

There are many situations in astrophysics, such as the bulk hypersonic 
motion of gas, where the kinetic energy of a fluid can dominate its internal 
energy by many orders of magnitude. In these situations, limitations on ma- 
chine precision can cause significant inaccuracy in the calculation of pressures 
and temperatures in the baryon gas. In order to address this issues, Enzo 
solves both the internal gas energy equation and the total energy equation 
everywhere on each grid, at all times. This dual energy formalism ensures 
that the method yields the correct entropy jump at strong shocks and also 
yields accurate pressures and temperatures in cosmological hypersonic flows. 

As a check on our primary hydrodynamic method, we also include an 
implementation of the hydro algorithm used in the Zeus astrophysical code. 

This staggered grid, finite difference method uses artificial viscosity as 
a shock-capturing technique and is formally first-order-accurate when using 
variable timesteps (as is common in structure formation simulations), and is 
not the preferred method in the Enzo code. 

2.4 Additional Physics Packages 

Several physics packages are implemented in addition to dark matter and 
adiabatic gas dynamics. The cooling and heating of gas is extremely important 
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in astrophysical situations. To this extent, two radiative cooling models and 
several uniform ultraviolet background models have been implemented in an 
easily extensible framework. 

The simpler of the two radiative cooling models assumes that all species 
in the baryonic gas are in equilibrium and calculates cooling rates directly 
from a cooling curve assuming Z = 0.3 Zm - The second routine, developed by 
Abel, Zhang, Anninos & Norman 0, Il2j . assumes that the gas has primor- 
dial abundances (ie, a gas which is composed of hydrogen and helium, and 
unpolluted by metals), and solves a reaction network of 28 equations which 
includes collisional and radiative processes for 9 seperate species (H, iJ + ,He, 
He + , He ++ , H~ , H 2 , and e~. In order to increase the speed of the cal- 
culation, this method takes the reactions with the shortest time scales (those 
involving H~ and i?2~) an d decouples them from the rest of the reaction net- 
work and imposes equilibrium concentrations, which is highly accurate for 
cosmological processes. See Anninos et al. 01 an d Abel et al. [llj for more 
information. 

The vast majority of the volume of the present-day universe is occupied by 
low-density gas which has been ionized by ultraviolet radiation from quasars, 
stars and other sources. This low density gas, collectively referred to as the 
"Lyman-a Forest" because it is primarily observed as a dense collection of 
absorption lines in spectra from distant quasars (highly luminous extragalactic 
objects), is useful because it can be used to determine several cosmological 
parameters and also as a tool for studying the formation and evolution of 
structure in the universe (see for more information). The spectrum of 
the ultraviolet radiation background plays an important part in determining 
the ionization properties of the Lyman-a forest, so it is very important to 
model this correctly. To this end, we have implemented several models for 
uniform ultraviolet background radiation based upon the models of Haardt & 
Madau 0. 

One of the most important processes when studying the formation and 
evolution of galaxies (and to a lesser extent, groups and clusters of galaxies 
and the gas surrounding them) is the formation and feedback of stars. We use 
a heuristic prescription similar to that of Cen & Ostriker to convert gas 
which is rapidly cooling and increasing in density into star "particles" which 
represent an ensemble of stars. These particles then evolve collisionlessly while 
returning metals and thermal energy back into the gas in which they formed 
via hot, metal-enriched winds. 

As mentioned in Section Enzo can be downloaded from the web at 
http : / / cosmos . ucsd . edu/enzo/ Vigorous code development is taking place, 
and we are in the process of adding ideal magnetohydrodynamics and a flux- 
limited radiation diffusion scheme to our AMR code, which will significantly 
enhance the capabilities of the code as a general-purpose astrophysical tool. 
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3 Adaptive Mesh Characteristics 

The adaptive nature of grid cells in the AMR simulations results in a wide 
range of baryon mass scales being resolved. Figure [21 shows the distribution 
of cells as a function of overdensity for a range of Enzo simulations in a sim- 
ulation volume which is 3 h~ l megaparsecs on a side. These simulations use 
either a 64 3 or 128 3 root grid and either 5 or 6 levels of refinement (such 
that Lbox/e = 4096, where Lbox is the box size and e is the smallest spatial 
scale that can be resolved). All grids are refined by a factor of 2.0, and grids 
are refined when dark matter density (baryon density) exceeds a factor of 4.0 
(2.0) times the mean density of cells at that level. In addition, a simulation is 
performed where the overdensity threshold is doubled. Initial conditions are 
generated using power spectra and methods common to cosmological simu- 
lations. Examination of Figure [3 shows that the entire density range in the 
simulations is covered by large numbers of cells. In particular, cells at low 
densities are well-resolved in these simulations, which is in stark contrast 
to simulations performed using Lagrangian methods, which are typically un- 
dersampled at low density. Raising the overdensity threshold for refinement 
decreases the total number of cells but their relative distribution as a function 
of overdensity is unchanged. In all simulations the total number of cells at the 
end of the run has increase by a factor of ~ 8 — 10 from the number of cells 
in the root grid. 

FigurcOHshows the distribution of number of cells at the end of a simulation 
as a function of the mass of baryons in that cell. Arrows indicate the mean cell 
mass contained on the root grid at the onset of the simulation for simulations 
covering the same spatial volume as the simulations described above with a 
64 3 , 128 3 or 256 3 root grid (labelled N64, N128 and N256, respectively). Over 
the course of the simulation the mean mass resolution, as indicated by the 
peak of the distribution, increases by almost an order of magnitude relative 
to the initial mass resolution, though the distribution of cell masses is quite 
large. Figure [3] shows that the mean cell mass as a function of overdensity 
(at the end of the simulation run) stays fairly constant, which lower mean cell 
masses in underdense regions and higher mean cell masses in highly overdense 
regions (presumably due to the limitation on the number of levels of adaptive 
mesh refinement allowed). The mean cell mass over the entire density range 
is between ^5 — 10 times better than the starting mass resolution for all sim- 
ulations. Runs with lower overdensity criteria for refinement have somewhat 
better mass resolution overall. 



4 Summary 

In this paper we have presented Enzo a cosmology code which combines colli- 
sionless N-body particle dynamics with a hydrodynamics package based on the 
piecewise parabolic method, all within a block-based adaptive mesh refinement 
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Fig. 2. Number of cells as a function of baryon overdensity (normalized by bin size) 
for a representative suite of Enzo cosmology simulations. The simulations are labelled 
such that the first number corresponds to the number of dark matter particles 
and the second number shows the root grid size, ie, 64dm/128grid means that the 
simulation has 64 3 dark matter particles and a 128 3 root grid. Simulations with two 
different hydrodynamic methods are used. PPM: Piecewise Parabolic Method. 
Zeus: The hydro method used in the Zeus astrophysical code. 



algorithm. Several other physics packages are implemented, including multi- 
ple models for gas cooling and ionization, a uniform ultraviolet background 
model for gas heating, and a prescription for star formation and feedback. 

Enzo is being released to the public as a community astrophysical simula- 
tion code. This code is being modified and documented to be as widely useful 
as possible, and can be found at http://cosmos.ucsd.edu/enzo/ Active 
development is taking place, centering on the addition of magnetohydrody- 
namics and a diffusive radiative transfer algorithm. 

Further information concerning the performance of the Enzo code (includ- 
ing a package of performance monitoring and visualization tools) is described 
in a contribution by James Bordner in this volume, and on the Enzo website. 
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Fig. 3. Number of cells as a function of baryonic mass (normalized to bin size) for 
several Enzo simulations with Lbox = 3k 1 Mpc. The arrows correspond to the mean 
mass resolution of the root grid for simulation volumes with the same volume but 
root grids with 64 3 , 128 3 or 256 3 cells (labelled N64, N128 and N256, respectively). 
See Figure |21 for a description of the line labels. 
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Fig. 4. Mean baryonic mass in cells as a function of baryon overdensity (normalized 
to bin size) in Enzo simulations. These are the same simulations (with the same 
labels) as in Figure^ Horizontal lines correspond to the initial mean mass resolution 



of simulations with the same volume but 64 3 , 
N64, N128 and N256, respectively). 



128 3 and 256 3 root grid cells (labelled 
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